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A  Random  Finite  Set  Approach  to  Space  Junk  Tracking  and 

Indentification 

Ba-Ngu  Vo\  Ba-Tuong  Vo\ 

^Department  of  Electrical  and  Computer  Engineering 
Curtin  University,  Perth,  WA,  Australia. 


ABSTRACT 

The  report  summarizes  the  feasibility  study  performed  for  AFRL/AOARD  regarding  the  use  of  Random  Finite 
Set  (RFS)  hltering  approach  for  tracking  orbital  space  debris.  Specifically,  we  investigate  the  capability  of  the 
RFS  approach  to  accommodate:  nonlinear  motion  and  measurement  models;  unknown  and  time  varying  target 
number;  and  limited  sensor  held-of-view.  This  report  shows  the  capability  of  the  random  hnite  set  approach 
to  provide  large  scale  multi-target  tracking.  In  particular  it  is  shown  that  an  approximate  filter  known  as  the 
labeled  multi-Bernoulli  filter  can  simultaneously  track  one  thousand  five  hundred  targets  in  clutter  on  a  standard 
laptop  computer. 


1.  INTRODUCTION 

Space  junk  is  the  term  commonly  used  to  refer  to  the  many  millions  of  pieces  of  waste  products  and  defunct 
objects  currently  in  orbit  around  the  earth,  which  are  largely  a  result  of  human  activities  in  relation  to  space 
exploration.^’®  Although  the  size  and  mass  of  individual  fragments  varies  from  millimetres  to  metres  and  grams 
to  tonnes,  even  the  most  seemingly  insignificant  pieces  can  pose  a  real  risk  of  collision  and  destruction  for  space 
related  infrastructure  and  exploration.  However  it  is  currently  understood  that  the  most  dangerous  debris  are 
those  with  diameters  between  1  and  10  centimetres,  of  which  it  is  estimated  that  there  are  several  hundred 
thousand,  since  they  are  small  enough  to  be  difficult  to  track  but  large  enough  to  cause  significant  damage. 
Fragments  larger  than  10cm  can  generally  be  tracked  and  avoided  for  time  being,  and  those  smaller  than  1cm 
can  generally  be  shielded  from  with  appropriate  apparatus. 

A  prerequisite  for  space  debris  management  is  to  understand  and  maintain  awareness  of  orbital  space  objects 
and  the  space  environment,  which,  in  the  space  domain,  is  called  Space  Situational  Awareness  (SSA).  Utilising 
available  data  from  multiple  sources  to  keep  track  of  multiple  space  objects  is  a  core  component  of  SSA.^  There 
are  three  major  approaches  to  tracking  multiple  objects:  Multiple  Hypotheses  Tracking  (MHT);^  Joint  Proba¬ 
bilistic  Data  Association  (JPDA);®  and  Random  Finite  Set  (RFS).^^  The  comprehensive  study  of  astrodynamics 
standards  in,^  identified  data  association,  information  fusion  and  nonlinear  estimation  as  three  key  research  areas 
for  tracking  multiple  space  objects.  MHT  and  its  variations  involve  the  propagation  of  data  association  hypothe¬ 
ses  in  time.  The  JPDA  approach  weights  individual  observations  by  their  association  probabilities.  MHT  and 
JPDA  are  classical  approaches  that  dominated  the  field  of  multi-target  tracking  and  are  well  documented  in.®’® 
The  RFS  approach  provides  a  general  systematic  treatment  of  multi-target  system  by  modeling  the  multi-target 
state  as  an  RFS  and  accommodates 

•  nonlinear  motion  and  measurement  models 

•  unknown  time  varying  target  number 

•  multi-sensor  data  with  support  for  heterogeneous  sensor  types 

•  limited  sensor  held-of-view 

Scalable  multi-target  tracking  solutions  are  needed  in  SSA  where  the  number  of  targets  is  large.®  However, 
multi-target  tracking  is  intrinsically  an  NP-hard  problem  and  the  complexity  of  multi-target  tracking  solutions 
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usually  do  not  scale  gracefully  with  problem  size.  The  RFS  approach  provides  the  flexibility  to  derive  multi-target 
tracking  algorithms  that  bypasse  the  data  association  problem  e.g  the  Probability  Hypothesis  Density  (PHD)/° 
Cardinalized  PHD/^  and  Multi-Bernoulli  filters. ^^^6  Although  RFS  approaches  have  been  extensively  applied  in 
aerospace  or  the  traditional  defence  related  problems,  there  has  been  very  little  activity  in  the  space  situational 
awareness  and  specifically  in  the  space  junk  monitoring  domain.  Current  works  include  an  application  to  tracking 
in  the  0-1  case,  or  Bernoulli  filtering,  where  at  most  one  target  is  present  with  signihcant  nonlinearities.'’ ’ 
Additionally  a  conceptual  solution  has  been  derived  assuming  a  fixed  a  known  number  of  targets.'^  Still  the 
feasibility  of  RFS/FISST  approaches  for  space  junk  monitoring  remains  an  open  question.  The  broad  objectives 
of  this  study  are  to 

•  Investigate  feasibility  of  RFS  based  approach  for  large  scale  real  time  space  junk  monitoring 

•  Investigate  methodological  and  computational  strategies  for  practical  implementation 

The  report  is  organized  as  follows.  Background  material  on  the  RFS  approach  to  multi-target  tracking  needed 
for  this  work  is  given  in  Section  2.  In  Section  3  we  present  the  labeled  multi-Bernoulli  filter-an  approximation 
to  the  optimal  Bayes  multi-target  tracking  filter-that  has  the  potential  to  address  large  scale  multi-target  prob¬ 
lems. In  Section  4  we  present  a  numerical  example  involving  1500  targets.  To  the  best  of  our  knowledge  this  is  the 
only  algorithm  that  can  track  such  a  large  number  of  targets  in  clutter  using  only  a  standard  laptop  computer. 

2.  LABELED  RANDOM  FINITE  SET  AND  MULTI-TARGET  TRACKING 

The  random  finite  set  (RFS)  or  finite  set  statistics  (FISST)  paradigm^^  for  multi-object  estimation  is  an  ideal 
candidate  platform  for  developing  top-down  algorithmic  solutions  for  space  junk  monitoring.  The  RFS  paradigm 
was  initially  developed  as  a  completely  new  paradigm  for  multi-object  estimation,  encompassing  a  principled 
mathematical  foundation  as  well  as  appropriate  calculus  like  tools  for  manipulating  probability  densities  of 
finite-set-valued  random  variables.  The  framework  is  general  enough  to  accommodate  non-linear  motion  and 
measurement  models,  including  image  measurement  models,  unknown  and  time  varying  number  of  target,  arbi¬ 
trary  sensor  field-of-view,  and  multi-sensor  data  support  for  heterogeneous  sensor  types  (and  even  non-standard 
data),  see.^^  Thus  the  RFS  approach  is  ideally  suited  to  networked  multi-sensor  multi-object  estimation  prob¬ 
lems,  such  as  space  situational  awareness. 

RFS-based  multi-target  filters  such  as  the  Probability  Hypothesis  Density  (PHD),^°  Cardinalized  PHD 
(CPHD)^^  and  multi-Bernoulli  filters^®’ can  avoid  data  association.  Indeed,  under  low  clutter,  it  was  demon- 
straeted  that  a  version  of  CPHD  filter  can  simultanenously  handle  over  a  thousand  targets  on  a  standard  laptop 
computer.  However  these  filters,  in  principle,  are  not  multi-target  trackers  because  they  rest  on  the  premise  that 
targets  are  indistinguishable. 

In,^”^  the  notion  of  labeled  RFSs  is  introduced  to  address  target  trajectories  and  their  uniqueness.  Moreover, 
an  analytic  solution  to  the  Bayes  multi-target  tracking  filter  known  as  the  d-generalized  labeled  multi-Bernoulli 
(d-GLMB)  hlter  was  developed.  Apart  from  producing  tracks,  the  d-GLMB  filter  does  not  assume  high  detection 
probability  and  low  clutter  and  outperforms  the  PHD/GPHD  and  multi-Bernoulli  hlters.  In  this  report  we  show 
that  the  labeled  RFS  technique  is  capable  of  tracking  simultaneously  thousands  of  targets  in  heavy  clutter  with 
different  birth  and  death  times,  on  a  standard  laptop  computer.  The  key  innovation  is  the  family  of  labeled  RFS 
conjugate  priors  in^^  and  a  principled  and  efficient  and  highly  parallelizable  approximation,  called  the  labeled 
multi-Bernoulli  filter.^® 

2.1  Random  Finite  Sets 

In  a  multiple  target  system,  the  number  of  targets  varies  with  time  due  to  the  appearance  and  disappearance 
of  targets,  and  the  number  of  measurements  received  at  each  time  step  does  not  necessarily  match  the  number 
of  targets  due  to  missed  detections  and  clutter.  The  objective  of  multiple  target  filtering  is  to  jointly  estimate 
the  number  of  targets  and  their  states  from  the  accumulated  observations.  Suppose  that  at  time  k,  there 
are  N{k)  target  states  Xk,i,  ■  ■  ■  ,Xk^N(k)j  each  taking  values  in  a  state  space  X  C  K"*,  and  M{k)  observations 
-Zfc.ij  •  ■  • )  Zk^M{k)  each  taking  values  in  an  observation  space  Z  C  Define  the  finite  sets 

Afc  ■  5  G  X, 

Xfc  =  •  7  G  Z, 
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to  be  the  multiple  target  state  and  multiple  target  observation  respectively. 

In  the  Bayesian  estimation  paradigm,  the  state  and  measurement  are  treated  as  realizations  of  random 
variables.  Since  the  (multi- object)  state  Xk  is  a  finite  set,  the  concept  of  a  random  finite  set  (RFS)  is  required 
to  cast  the  multi-object  estimation  problem  in  the  Bayesian  framework.  An  RFS  is  simply  a  random  variable 
that  take  values  as  (unordered)  finite  sets,  i.e.  a  finite-set-valued  random  variable.  The  space  of  finite  subsets 
of  X  does  not  inherit  the  usual  Euclidean  notion  of  integration  and  density.  Hence,  standard  tools  for  random 
vectors  are  not  appropriate  for  RFSs.  Mahler’s  Finite  Set  Statistics  (FISST)  provides  powerful  yet  practical 
mathematical  tools  for  dealing  with  RFSs.^°A2  notion  of  probability  densities,  generating  functionals,  and 
calculus  like  tools  for  are  provided  by  the  theory  of  finite  set  statistics.^^’^'^  For  connections  between  aspects  of 
the  FISST  notion  of  probability  densities  and  standard  measure  theoretic  probability  theory  we  refer  the  reader 
t0.14 


2.2  Labeled  Random  Finite  Sets 

To  incorporate  target  identity,  each  state  a;  G  X  is  augmented  with  a  unique  label  £  G  L=  {£i  :  i  £  N},  where  N 
denotes  the  set  of  positive  integers  and  the  £i^s  are  distinct.  More  detail  on  labeled  RFS  can  be  found  in.^^’^® 

Definition  2.1.  A  labeled  RFS  with  state  space  X  and  (discrete)  label  space  L  is  an  RFS  on  XxL  such  that 
each  realization  has  distinct  labels. 

Let  C  :  XxL  — >■  L  be  the  projection  C{{x,  £))  =  £,  then  a  finite  subset  set  X  of  XxL  has  distinct  labels  if  and 
only  if  X  and  its  labels  £(X)  =  {£(x)  :xGX}  have  the  same  cardinality,  i.e.  (I|x| (|>C(X)|)  =  1.  The  function 
A(X)  =  (5|x| (|>C(X)|)  is  called  the  distinct  label  indicator. 

The  unlabeled  version  of  a  labeled  RFS  is  obtained  by  simply  discarding  the  labels.  Consequently,  the 
cardinality  distribution  (the  distribution  of  the  number  of  objects)  of  a  labeled  RFS  the  same  as  its  unlabeled 
version. 

For  the  rest  of  the  paper,  single-object  states  are  represented  by  lowercase  letters,  e.g.  x,  x  while  multi-object 
states  are  represented  by  uppercase  letters,  e.g.  X,  X,  symbols  for  labeled  states  and  their  distributions  are 
bolded  to  distinguish  them  from  unlabeled  ones,  e.g.  x,  X,  tt,  etc.,  spaces  are  represented  by  blackboard  bold 
e.g.  X,  Z,  L,  N,  etc.,  and  the  class  of  finite  subsets  of  a  space  X  is  denoted  by  J'(X).  The  inner  product  of  two 
functions  /  and  g  is  denoted  by  (/,(/)  =  /  f{x)g{x)dx.  The  integral  of  a  function  f  XxL  — >•  K  is  given  by 

f  f(x)dx  =  ^  f  t{{x,£))dx. 


2.3  Labeled  multi-Bernoulli  RFS 

A  labeled  multi-Bernoulli  RFS  X  with  state  space  X,  label  space  L  and  parameter  set  :  C  G  'Ll, 

is  a  multi-Bernoulli  RFS  on  X  augmented  with  labels  corresponding  to  the  successful  (non-empty)  Bernoulli 
components,  i.e.  if  the  Bernoulli  component  {r^‘’\p^‘='^)  yields  a  non-empty  set,  then  the  label  of  the  corresponding 
state  is  given  by  q;(C),  where  a  :  T  — >•  L  is  a  1-1  mapping.  A  labeled  multi-Bernoulli  RFS  has  (multi-target) 
density  given  by^^ 


7r({(xi,4), 


iXn,£n)})=6r.i\{£u-£n}\)  J]  H 

C6'p  i=i 


I  —  7*(^  ^  (D) 


where 


1,  if  A  C  y 
0,  otherwise. 


is  the  inclusion  function.  For  convenience  we  use  the  abbreviation  tt  =  for  the  density  of  a 

labeled  multi-Bernoulli  RFS.  Although  the  formulation  allows  for  a  general  mapping  a  for  the  labels,  in  this 
work  we  assume  that  a  is  an  identity  mapping  in  order  to  simplify  notations.  Thus  superscripts  for  component 
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indices  correspond  directly  to  the  label  in  question.  The  density  of  a  labeled  multi-Bernoulli  RFS  with  parameter 
set  TT  =  is  given  more  compactly  by 


where  =  nx£xP(x)> 


7r(X)  =  A(X)w(£(X))p^ 


w{L) 

pix,£) 


n(i-.«)n 

iSL  ieL 


1  -  r(^)  ’ 


p‘'^\x). 


(1) 

(2) 

(3) 


2.4  Generalized  Labeled  Multi-Bernoulli  RFS 

A  generalized  labeled  multi-Bernoulli  RFS  is  a  labeled  RFS  on  XxL  distributed  according  to 


7r(X)  =  A(X)^u;(^)(£(X)) 


cGC 


(c) 


where  C  is  a  discrete  index  set,  (L)  and  satisfy 


LCLcGC 

/p<=>(x,f)d. 


1, 

1. 


(4) 

(5) 

(6) 


The  cardinality  distribution  of  a  generalized  labeled  multi-Bernoulli  RFS  is  given  by 

LgJ='„(L)  cGC 


(7) 


The  labeled  multi-Bernoulli  RFS  distributed  by  (1)  is  a  special  case  of  the  generalized  labeled  multi-Bernoulli 
RFS  with 


p^‘^\x,i)  =  p^^\x) 

iGL-L  ^GL 

comprising  a  single  component  in  which  case  the  superscript  (c)  is  omitted. 


2.5  GLMB  Multi-target  Tracking  Filter 

In  labeled  RFS  multi-target  tracking,  targets  are  identified  by  an  ordered  pair  of  integers  i  =  {k,i),  where  k 
is  the  time  of  birth  and  z  €  N,  with  N  denoting  the  set  of  natural  numbers,  is  a  unique  index  to  distinguish 
objects  born  at  the  same  time.  The  label  space  for  objects  born  at  time  k,  denoted  as  L^,  is  then  {k}  x  N.  An 
object  born  at  time  k,  has  state  x  G  XxL^.  The  label  space  for  targets  at  time  k  (including  those  born  prior  to 
A:),  denoted  as  Lo:fc,  is  constructed  recursively  by  Lqiz,  =  Lo:fe_i  U  (note  that  Lo:fc_i  and  are  disjoint).  A 
multi-object  state  X  at  time  k,  is  a  finite  subset  of  XxLorfc. 

Let  7rfc(-|Xfc)  denotes  the  multi-target  posterior  density  at  time  k,  and  denotes  the  multi-target  pre¬ 

diction  density  to  time  k.  Then,  the  multi-target  Bayes  recursion  propagates  ttz,  in  time^^  according  to  the 
following  update  and  prediction 


k{^k  l-^fc)  — 


k\k—l{^k^ 


75fc(zfc|X)7rfc|fc_i(X)5X’ 


(8) 

(9) 
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where  gk{‘\‘)  is  the  multi-target  likelihood  function  at  time  k,  ffc+i|fe(-|-)  is  the  multi-target  transition  density,  to 
time  k-\-l,  and  the  integral  is  a  set  integral  defined  for  any  function  f  :  J^(XxL)  — >■  K  by 

J  f(X)(5X  =  ^^  J  f({xi,...,xi})d(xi,...,xi). 

i=0 

The  multi-target  posterior  density  captures  all  information  on  target  number,  and  individual  target  states. 
The  multi-target  likelihood  function  encapsulates  the  underlying  models  for  detections  and  false  alarms  while 
the  multi-target  transition  density  encapsulates  the  underlying  models  of  target  motions,  births  and  deaths. 

For  convenience,  in  what  follows  we  omit  explicit  references  to  the  time  index  fc,  and  denote  L  =  Lo:fc, 
B  =  Lfe_|_i,  L_|_=  L  U  B,  Tv~TVk,  7r+=:7rfc_|_i|fc,  g=gk,  In  this  work  we  assume  the  space  of  birth  label  B  is 

hnite. 


2.6  GLMB  Prediction 

Given  the  current  multi-object  state  X,  each  state  {x,  £)  G  X  either  continues  to  exist  at  the  next  time  step  with 
probability  ps{x,£)  and  evolves  to  a  new  state  (x+,£+)  with  probability  density  f{x+\x,£)6i{£+),  or  dies  with 
probability  1  —  ps{x,£).  The  set  of  new  targets  born  at  the  next  time  step  is  distributed  according  to 

fB(Y)  =  A(YVB(/:(Y))[pBr  (10) 

The  birth  density  is  defined  on  X  x  B  and  fs  (Y)  =  0  if  Y  contains  any  element  y  with  £(y)  ^  B.  The  birth 
model  (10)  covers  both  labeled  Poisson  and  labeled  multi-Bernoulli.  For  a  labeled  multi-Bernoulli  birth  model: 

wsiJ)  =  n  (i“T)  (^^) 

iGB-J  i&J 

Pb{x,£)  =  Pb\^)-  (12) 

The  multi-target  state  at  the  next  time  X+  is  the  superposition  of  surviving  targets  and  new  born  targets. 
Assuming  that  targets  evolve  independently  of  each  other  and  that  births  are  independent  of  surviving  targets, 
it  was  shown  in^^  that  the  multi-target  transition  kernel  is  given  by 

f  (X+|X)  =  fs(X+  n  (X  X  L)|X)fB(X+  -  (X  X  L))  (13) 


where 


fs(W|X) 

d>(W;a;,£) 


A(W)A(X)l£(x)(/:(W))  [$(W;  •)]'' 

(  ps{x,£)f{x+\x,£),  if(x+,.^)GW 

\  l-ps{x,£),  ii£iCiW) 


Proposition  2.2.  If  the  current  multi- object  prior  is  a  generalized  labeled  multi- Bernoulli  of  the  form  (j),  then 
the  predicted  multi-object  density  is  also  a  generalized  labeled  multi- Bernoulli  given  by 


7r+(X+)  =  A(X^ 


'E' 

cGC 


<)(£(X+)) 


P+ 


(14) 


where 


w+\l) 

P+\x,£) 

p''s\x,£) 


WB{L-l,)w‘'g\Lnl,), 
^hi£)p^s^ix,£)  -k  (1  -  lTLi£))pB{x,£), 
{psi-,£)fix\-,£),p^^'i  {■,£)) 
pf{£) 


(15) 

(16) 

(17) 
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(18) 

(19) 


/CL 

4=)(£)  =  l-(^ps{;£),p(‘^\;i)).  (20) 


2.7  GLMB  Update 

For  a  given  multi-target  state  X,  at  time  k,  each  state  (x,£)  €  X  is  either  detected  with  probability  {x,£)  and 
generates  a  point  z  with  likelihood  g{z\x,£),  or  missed  with  probability  1— p£i(x,£).  The  multi-object  observation 
Z  =  {zi, ...,  z\z\}  is  the  superposition  of  the  detected  points  and  Poisson  clutter  with  intensity  function  k. 

Definition  2.3.  An  association  map  (for  the  current  time)  is  a  function  0  :  L  — >■  {0, 1, ...,  |Z|}  such  that 
6{{)  =  6{i')  >  0  implies  i  =  i! .  The  set  0  of  all  such  association  maps  is  called  the  association  map  space.  The 
subset  of  Q  with  domain  I  is  denoted  by  0(1). 

An  association  map  describes  which  tracks  generated  which  measurements,  i.e.  track  £  generates  measurement 
zgs^i)  £  Z,  with  undetected  tracks  assigned  to  0.  The  condition  d(i)  =  d(i')  >  0  implies  i  =  i',  means  that  a  track 
can  generate  at  most  one  measurement  at  any  point  in  time. 

Assuming  that,  conditional  on  X,  detections  are  independent,  and  that  clutter  is  independent  of  the  detec¬ 
tions,  the  multi-object  likelihood  is  given  by 


g(Z|X)  =  ^ 

eee(£(x)) 


where 


i>z{x,£;e)  = 


if0(£)>O 
1 -p_D(a;,.^),  if  6l(£)  =  0 


(21) 


(22) 


Proposition  2.4.  If  the  prior  distribution  is  a  generalized  labeled  multi- Bernoulli  of  the  form  ()),  then,  under 
the  multi-object  likelihood  function  (21),  the  posterior  distribution  is  also  a  generalized  multi- Bernoulli  given  by 


7r(X|Z)  =  A(X)^  ^ 


f’^^(£(X))  [p("’^)(-|Z) 


X 


cGC^G0(£(X)) 


(23) 


where 


p^^’^\x,£\Z) 


(c,e) 

Vz 


(^) 


oc 


p^^\x,£)il}z{x,£-,9) 

(p(^)(-,£),V'z(-,^;0)), 


(24) 

(25) 

(26) 


For  a  valid  label  set  L,  the  updated  weight  is  proportional  to  the  prior  weight  scaled  by  the 

product  of  single- object  normalizing  constants. 

The  above  results  are  the  prediction  and  update  step  of  the  (5-generalized  labeled  multi-Bernoulli  tracking 
filter, also  known  as  the  Vo- Vo  filter  (see  Mahlerbook2014).  Tractable  techniques  for  truncating  the  posterior 
and  prediction  densities  were  proposed  based  on  the  fc-shortest  paths  and  ranked  assignment  algorithms,  see^® 
for  further  details. 
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3.  LABELED  MULTI-BERNOULLI  FILTER  FOR  LARGE  SCALE  TRACKING 


The  labeled  multi-Bernoulli  filter  is  an  approximation  of  the  (5-generalized  labeled  multi-Bernoulli  filter  using 
labeled  multi-Bernoulli  RFS.  The  number  of  components  in  the  the  Vo- Vo  hlter  grows  exponentially,  whereas  the 
growth  is  linear  for  a  labeled  multi-Bernoulli  representation.  In  this  section  we  outline  the  labeled  multi-Bernoulli 
hlter  and  its  implementation,  further  detail  can  be  found  in.^® 


3.1  LMB  Prediction 

It  was  shown  in^®  that  the  labeled  multi-Bernoulli  is  closed  under  the  prediction  step. 


Proposition  3.1.  Suppose  that  the  multi-target  posterior  density  is  a  labeled  multi- Bernoulli  with  parameter 
set  TT  =  o,nd  that  the  multi-target  birth  model  is  a  labeled  multi-Bernoulli  with  parameter  set 

tvb  =  {rg\pB^}eeM  (where  LflB  =  0^,  then  the  multi-target  predicted  density  is  also  a  labeled  multi- Bernoulli 
with  parameter  set 


where 


-+  =  { 

(27) 

(28) 

P-i-.S 

=  {Ps{-,T)f{x\-,e),p{-,£))  /igs{£), 

(29) 

Vsi^) 

=  {ps{-,£),p{-,£)) 

(30) 

Remark  3.2.  The  multi-target  prediction  for  a  labeled  multi- Bernoulli  actually  coincides  with  the  prediction  for 
the  unlabeled  multi- Bernoulli  and  interpreting  the  component  indices  as  track  labels. 


3.2  LMB  Update 

While  the  labeled  multi-Bernoulli  family  is  closed  under  the  prediction  operation,  it  is  not  closed  under  the 
update  operation.  Inspired  by  the  multi-Bernoulli  hlter, we  seek  the  labeled  multi-Bernoulli  approximation  that 
matches  the  hrst  moment  of  the  multi-target  posterior  density.  One  such  labeled  multi-Bernoulli  approximation 

•  •  •  1 Q 

IS  given  m. 

Proposition  3.3.  Suppose  that  the  multi-target  predicted  density  is  a  labeled  multi- Bernoulli  with  parameter 
set  •  The  labeled  multi- Bernoulli  that  matches  exactly  the  first  moment  of  the  multi-target 

posterior  density  is  7r(-|Z)  =  where 


(i+,9)ej^(i.+)x0(i+) 


1 


w^^+'^\z)h4e)p^<^\x,e\z), 


(/+.0)GJ"(L+)xe(/+) 


w^^+’^\Z)  oc 


p^'^\x,i\Z)  = 


W 

Vz 


'^-1-  > 


iei+ 


p+{x,£)ipzix,T,6) 


(31) 

(32) 

(33) 

(34) 

(35) 


The  key  advantage  of  the  labeled  multi-Bernoulli  update  is  that  it  only  involves  one  approximation  of  the 
multi-target  posterior  density.  In  contrast  the  multi-Bernoulli  hlter  requires  two  approximations  on  the  multi¬ 
target  posterior  probability  generating  functional. 
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The  mean  cardinality  of  the  labeled  multi-Bernoulli  approximation  is  identical  to  that  of  the  full  posterior. 
However,  the  cardinality  distributions  themselves  differ,  since  the  cardinality  distribution  of  the  labeled  multi- 
Bernoulli  approximation  follows  the  one  for  a  multi  Bernoulli  RFS  while  the  cardinality  distribution  of  the 
full  posterior  follows  that  of  a  Generalized  labeled  multi-Bernoulli  (GLMB)  RFS.  Although  there  are  many 
possible  choices  of  approximate  labeled  multi-Bernoulli  posteriors,  the  above  choice  preserves  the  estimated 
spatial  densities  of  each  track  and  matches  exactly  the  first  moment.^® 

3.3  Efficient  Implementation  of  the  Labeled  Mnlti-Bernoulli  Filter 

The  labeled  multi-Bernoulli  form  can  naturally  be  decomposed  into  products  of  smaller  labeled  multi-Bernoullis 
that  allows  for  parallel  update  of  ‘groups’  of  closely  spaced  targets  and  their  associated  measurements,  which 
drastically  reduces  computation  time.  Target  grouping  is  based  on  a  standard  gating  procedure  which  also  parti¬ 
tions  the  measurement  set,  resulting  in  groups  of  targets  and  measurements  which  can  be  considered  statistically 
independent.®  Each  resultant  group  is  usually  much  smaller  than  the  entire  multi-target  state  and  measurement, 
and  can  then  be  updated  in  parallel,  which  is  usually  much  simpler  and  faster  than  updating  the  entire  multi¬ 
target  state  as  a  single  group.  The  update  for  each  group  is  performed  by  expanding  the  labeled  multi-Bernoulli 
prediction  into  d-GLMB  form,  and  performing  a  standard  d-GLMB  update  resulting  in  a  d-GLMB  posterior.®^ 
The  posterior  for  each  group  is  then  collapsed  to  a  matching  labeled  multi-Bernonlli  approximation  after  which 
groups  are  recombined  and  the  recursion  continues.  This  proposed  implementation  applies  to  both  Gaussian 
Mixture  and  Sequential  Monte  Carlo  based  representations  for  the  underlying  single  target  densities.^® 

3.3.1  Grouping  and  Gating 

Let  ...,L^^}  be  a  partition  of  the  label  set  L+  =  L  U  B,  and  be  a  partition  of  the 

measurement  set  Z.  A  grouping  is  defined  as  the  set  of  pairs  {(L^\  ...,  and  each  pair 

=  (L^\  Z*^”))  is  referred  to  as  a  group.  Note  that  Z*^°^  is  the  set  of  measurements  which  are  not  assigned 
to  any  target  labels. 

We  are  interested  in  groupings  of  target  labels  and  likely  corresponding  measurements.  To  generate  such 
groupings,  we  start  with  an  initial  grouping  of  each  labeled  Bernoulli  component  £  and  any  associated  measure¬ 
ments  which  fall  within  its  prediction  gates: 

=  m,  {z  :  z)  <  V7}).  (36) 

where  z)  is  the  Mahalanobis  distance  (MHD)  between  the  predicted  measurement  for  track  £  and  the 

received  measurement  z  &  Z,  and  7  is  the  gating  distance  threshold  calculated  using  the  inverse  Chi-squared 
cumulative  distribution  corresponding  to  the  desired  cr-gate  size  for  gating  of  measurements  from  tracks.  Then, 
tentative  groups  with  common  measurements 

z(d  n  z(^')  ^  0  (37) 

are  merged  as  follows 

=  (L^  UL^'^ZW  U  Z(^')).  (38) 

The  merging  is  repeated  for  all  tentative  groups  until  there  are  no  common  measurements.  Finally,  a  total  of  N 
groups  g^^\  ...,  g^^'>  of  tracks  and  associated  measurements  is  obtained.  Consequently,  the  predicted  multi-target 
density  can  be  rewritten  as 

N 

7^+  =  U  (39) 

where 


Distribution  A:  Approved  for  pubiic  reiease.  Distribution  is  Uniimited. 


3.3.2  Parallel  5-GLMB  Group  Updates 


Since  the  predicted  multi-target  density  for  each  group  is  a  labeled  multi-Bernoulli  of  the  form  (39),  we  need 
to  express  it  in  (5-GLMB  form,  in  order  to  apply  the  data  update.  The  J-GLMB  form  for  the  f-th  group 
g(i)  =  is  given  by 


7r«(X+)  =  A(X+)  ^  u;;^+)<5,^(£(X+))b+]^+  (41) 


where 


w 


(1+) 

-\-4 


n  (i-4‘’)  n 


Due  to  the  smaller  label  space  within  one  group,  the  number  of  hypotheses  across  all  groups  is  significantly 
smaller  than  for  the  number  of  hypotheses  in  case  of  a  single  big  group. 

A  brute-force  way  to  enumerate  the  sum  is  to  generate  all  possible  combinations  for  a  set  of  labels 
and  cardinalities  n  =  0, 1, . . . ,  |L^^|.  The  number  of  combinations  for  each  cardinality  is  given  by  the  binomial 
coefficient  G(|L^^|,n)  =  |L^^|!/(n!(|L^^|  — n)!)  and  the  number  of  combinations  for  a  set  of  track  labels  is 

Consequently  explicit  evaluation  of  all  combinations  is  only  possible  for  small  For  larger  |L^|!^|,  the  sum  can 

be  approximated  to  its  k  most  significant  terms  by  use  of  the  /c-shortest  paths  algorithm  without  enumerating 
all  possible  terms. Consequently,  1+  only  consists  of  the  most  significant  hypotheses. 

The  (5-CLMB  update  for  each  group  i  is  given  by 


rW 


(X+|ZW)=A(X+) 


E 


ri;(^+’^)(Z(*))5/^(£(X+))  [pW(-|zW) 


(/+.e)GJG(L«)xej+ 


where 


n 


p+{x,l)'ipz(i'i  0) 


n  (i-4'>)  n  c-> 


PD{x,e)pGg{zi'(\~,\x, 


if6i(£)>0 
l-pDix,£)pG,  if  6»(£)  =  0 


(42) 

(43) 

(44) 

(45) 


Observe,  that  (45)  has  to  incorporate  the  gating  probability  pc,^  since  small  gates  increase  the  probability  of 
missed  detection. 

By  representing  the  complete  predicted  distribution  using  several  groups,  the  number  of  track  labels  | 
within  each  group  is  significantly  lower  than  the  total  number  of  possible  labels  |L+|.  Since  the  update  is 
combinatorial,  the  number  of  components  or  hypotheses  within  each  group  grows  exponentially  in  the  number 
of  track  labels  |L^^|.  Thus,  for  large  |L^^|  a  truncation  of  the  i5-CLMB  distribution  is  required.  The  truncation 
can  be  realized  using  ranked  assignment  algorithm  which  evaluates  only  the  M  most  significant  hypotheses 
without  evaluating  all  possible  solutions.^®  Since  the  complexity  of  ranked  assignment  algorithm  is  cubic,  the 
computational  costs  for  the  evaluation  of  multiple  groups  is  generally  cheaper,  compared  to  the  evaluation  for  a 
single  large  group.  Moreover,  the  evaluation  for  each  group  can  be  performed  in  parallel. 
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3.3.3  Approximation  of  the  updated  (5-GLMB  as  labeled  multi-Bernoulli 

After  performing  the  group  updates  across  all  groups  for  i  =  the  (5-GLMB  form  for  each  group  is 

then  collapsed  back  to  labeled  multi-Bernoulli  form 

7rW(-|Z«)  « 

where  are  calculated  according  to  (31)-(32)  and  taking  the  union  of  the  approximate  labeled  multi- 

Bernoulli  groups  given  by  yields  the  labeled  multi-Bernoulli  approximation  to  the  multi-target  posterior 

N 

7r(.|Z)«5f(.|Z)  =  U5f«(-|^^'^)-  (46) 

i=l 

Remark  3.4.  Theoretically  in  the  worst  case  where  all  targets  are  close  together,  it  may  not  be  possible  to 
partition  into  smaller  groups.  In  this  case,  the  complexity  of  the  update  step  is  the  same  as  that  of  the  S-GLMB 
filter.  However,  our  experience  with  empirical  data  suggests  that  in  most  scenarios  we  can  partition  the  targets 
into  many  groups  with  small  numbers  of  targets. 

3.4  Track  Extraction 

Since  tracks  are  represented  after  update  by  a  labeled  multi-Bernoulli,  an  obvious  track  extraction  scheme  is  to 
pick  all  tracks  with  an  existence  probability  higher  than  an  application  specific  threshold: 

X=  >r?}  (47) 

where  x  =  arg^,  maxp(^^(x).  On  the  one  hand,  choosing  a  high  value  for  d  significantly  reduces  the  number 
of  clutter  tracks  at  the  cost  of  a  delayed  inclusion  of  a  new-born  track.  On  the  other  hand,  low  values  for  D 
report  new-born  tracks  immediately  at  the  cost  of  a  higher  number  of  clutter  tracks.  Choosing  a  high  value  for 
D,  additional  care  has  to  be  taken  for  missed  detections.  In  case  of  pu  «  1,  a  missed  detection  considerably 
reduces  the  existence  probability.  Consequently,  one  missed  detection  might  suppress  the  output  of  a  previously 
confirmed  track  with  m  1. 

To  mitigate  this  issue,  a  hysteresis  is  used,  which  activates  the  output  if  the  maximum  existence  probability 
Umax  of  a  track  i  has  once  exceeded  an  upper  threshold  du  and  if  the  current  existence  probability  is  higher 
than  a  lower  threshold  di: 

X  =  {(x,£):r^L>^uAr(^)  (48) 

4.  A  LARGE  SCALE  TRACKING  DEMONSTRATION 

This  demonstration  of  the  labeled  multi-Bernoulli  (LMB)  filter  involves  3D  states  and  observations  in  x,  y,  z 
coordinates,  and  features  1500  targets  and  various  births  and  deaths  throughout.  Figure  1  shows  the  true  tracks 
in  3D-space  to  indicate  the  density  and  number  of  targets  in  this  scenario.  We  note  the  the  implementation  for 
this  example  is  a  standard  serial  implementation  and  has  not  been  parallelised.  Parallelisation  would  significant 
improve  teh  speed  of  the  as  well  as  scalability. 

The  following  dynamical  and  measurement  models  are  used.  The  target  state  variable  is  a  vector  of  3D  position 
and  velocity  Xk  =  [  Px,k,Px,k,Py,k,Py,k,Pz,k,Pz,k  ]^-  The  single-target  transition  model  is  linear  Gaussian  with 
transition  density  /fc|fc-i(a;|C)  =  J^(x]Fk-iC,,Qk-i)  where 
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A  =  Is  is  the  sampling  period,  and  Ui,  =  5km/ is  the  standard  deviation  of  the  process  noise.  The  probability 
of  survival  ps,k  =  0.99.  The  cardinality  distributions  involved  have  a  maximum  support  of  Amax  =  2000  terms. 
The  birth  process  is  modeled  as  a  Poisson  RFS  with  an  average  of  15  births  per  scan,  coming  from  3  different 
positions  in  the  3D  surveillance  region. 

The  probability  of  detection  pD,k  =  0.98.  The  single-target  measurement  model  is  also  linear  Gaussian  with 
likelihood  gk{z\x)  =  HkX,  Rk)  where 
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(Je  =  10km  is  the  standard  deviation  of  the  measurement  noise.  Clutter  is  modeled  as  a  Poisson  RFS  with 
uniform  distribution  at  an  average  of  100  points  per  scan. 


Figure  1.  3D  tracks  for  approximately  1500  targets. 


Note  that  we  use  the  linear  Gaussian  model  for  this  example  as  a  proof-of-concept,  to  show  that  the  proposed 
LMB  filter  can  handle  a  large  number  of  targets.  Non-linear  models  that  include  non-uniform  clutter,  non- 
uniform  sensor  field-of-view  can  be  easily  accommodated  using  a  particle  implementation.  However,  the  cost 
of  single-target  filtering  increases,  which  in  turn  increase  the  computational  requirement  for  the  multi-target 
tracker. 


For  performance  evaluation,  the  Optimal  Sub-Pattern  Assignment  (OSPA)  distance  between  the  estimated 
and  true  multiple  target  state  is  used  as  the  estimation  error. The  OSPA  metric  is  defined  as  follows.  Let 
d/'^\x,  y)  :=  min  (c,  ||a;  —  y\\)  for  x,y  G  X,  and  H^  denote  the  set  of  permutations  on  {1,  2, . . . ,  fc}  for  any  positive 
integer  k.  Then,  for  p  >  1,  c  >  0,  and  X  =  {xi, . . . ,  Xm}  and  Y  =  {yi, . . . ,  p„}, 


4^)(A,y)  := 


min 


^  d^‘'\xi,  -k  cP{n-m) 


1 

V 


if  TO  <  n,  and  d!'p\x,Y)  :=  d^\Y,X)  if  to  >  n;  and  d^\x,Y)  =  (^\y,X)  =  0  if  to  =  n  =  0.  The  OSPA 
distance  is  interpreted  as  a  p-th  order  per-target  error,  comprised  of  a  p-th  order  per-target  localization  error 
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and  a  p-th  order  per-target  cardinality  error.  The  order  parameter  p  determines  the  sensitivity  of  the  metric  to 
outliers,  and  the  cut-off  parameter  c  determines  the  relative  weighting  of  the  penalties  assigned  to  cardinality 
and  localization  errors. 

Plots  of  true  and  estimated  multiple  target  cardinality  as  well  as  the  OSPA  miss-distance  (for  p  =  1  and 
c  =  50km)  from  the  labeled  multi-Bernoulli  filter  are  shown  in  figures  2  and  3.  From  the  plot  of  the  cardinality  in 
hgure  2,  it  appears  that  the  labeled  multi-Bernoulli  filter  determines  the  number  of  targets  within  a  reasonable 
error.  The  plots  of  the  OSPA  miss  distance, from  figure  3,  confirm  the  correct  operation  of  the  filter;  the  errors 
are  consistent  with  the  measurement  noise  in  the  model  and  the  number  of  target  as  seen  from  the  cardinality 
plot.  Execution  time  per  cluster  for  each  iteration  of  the  filter  on  a  single  standard  laptop  computer  is  shown  in 
figure  4.  Consequently,  with  modest  computational  resources,  the  labeled  multi-Bernoulli  filter  can  potentially 
deliver  real  time  performance  in  large  scale  tracking. 


Figure  2.  LMB  estimated  cardinality  versus  time 


Time 

Figure  3.  OSPA  miss  distance  given  by  total,  localization  component  and  cardinality  component 


5.  CONCLUSION 

The  RFS  approach  is  ideally  suited  to  networked  multi-sensor  multi-object  estimation  problems,  such  as  space 
situational  awareness.  It  provides  a  general  systematic  treatment  of  multi-target  system  and  accommodates: 
nonlinear  motion  and  measurement  models;  unknown  time  varying  target  number;  multi-sensor  data  with  support 
for  heterogeneous  sensor  types;  limited  sensor  held-of-view.  More  importantly,  it  provides  scalable  algorithms 
that  can  be  uses  to  track  a  large  number  of  targets.  A  3D  scenario  using  Gaussian  mixture  implementations  has 
shown  that  the  labeled  multi-Bernoulli  filter  can  simultaneously  track  up  to  1500  targets,  with  different  birth 
and  death  times,  on  a  single  standard  laptop  computer.  Performance  analysis  of  the  hlter  also  suggests  that  it 
accurately  estimates  target  numbers  and  states  in  dense  situations.  We  note  the  the  implementation  for  this 
example  is  a  standard  serial  implementation  and  has  not  been  parallelised.  Without  parallelization,  execution 
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Figure  4.  Average  execution  time  per  cluster  for  each  iteration 

times  for  each  iteration  was  of  the  order  of  tens  of  seconds.  Parallelisation  would  significantly  improve  the 

speed  of  the  algorithm  as  well  as  scalability.  Consequently,  with  modest  computational  resources,  the  labeled 

multi-Bernoulli  filter  can  potentially  scale  to  real  time  large  scale  tracking. 
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